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ABSTRACT 

In this paper a new efficient algorithm for computation of radio wave ray trajectories in isotropic 
media is described. The algorithm is based on an original second-order difference scheme with a specific 
"length-conservation" property, which allows to resolve the ray shape even in the regions where its 
curvature is high. Besides the scheme, the algorithm includes a number of mechanisms securing its 
correctness and stability. The algorithm is intended for obtaining multi-pixel solar images in wide range 
of radio frequencies, and it is designed to be used in studies of the solar environment with modern 
high-resolution radiointerferometers and radiotelescopes such as the Murchison Widefield Array (MWA). 



1. Introduction 

The problems of image simulation and recon- 
struction in solar radioastronomy require compu- 
tation of multiple ray trajectories. A typical image 
of the sun consists of thousands to millions pixels, 
and the brightness of each of the pixels is obtained 
as an integral over an individual ray trajectory. 
Due to non-uniform distribution of the index of 
refraction in heliospheric plasmas at lower frequen- 
cies (up to a few GHz) the ray trajectories cannot 
be treated as straight lines. The ray-tracing prob- 
lems require significant amount of computations, 
and, hence, time. It is possible to considerably 
reduce the computation time through the use of 
efficient ray-tracing algorithms. 



In this work a ray tracing algorithm is pre- 
sented. It allows one to calculate trajectories of 
electromagnetic rays in space plasmas with speci- 
fied density distributions. In most applications the 
radio rays are traced to calculate certain physical 
quantities, such as intensity or black body bright- 
ness temperature, through path integration over 
the ray trajectories. In its software implementa- 
tion the algorithm performs the integration and 
the ray trajectory calculation simultaneously. 

We start with a brief review of the facts from 
the physics of electromagnetic wave propagation 
in plasmas pertinent to the ray tracing. We show 
the relationship between the plasma density, the 
dielectric permittivity, and the index of refraction 
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as well as between their gradients. In the next sec- 
tion we derive the ray equations from the Fermat's 
principle. We transform the set of ray equations 
so that their solution is a naturally parameterized 
curve. This representation reveals the mathemat- 
ical analogy of the ray equations and the equa- 
tions for a particle motion in magnetic field. Us- 
i ng th is analogy, we derive a Crank fc NicolsonI 
(1 19471 ) secon d-order difference scheme using the 
Boris (Il970l) approach. In the Boris' scheme the 
particle energy, or the squared velocity, is con- 
served. For the naturally parameterized curve like 
that we use here the first derivative of the vector 
to be solved from the equation is not velocity (as 
present in the equation of the particle motion), 
but the unity-length direction cosine vector. Our 
scheme conserves the squared direction vector at 
unity. To achieve this, we employ the ideas which 
are often used while simulating the particle mo- 
tion and which are known to conserve the velocity 
vector length in the magnetic field. 

The difference scheme coding is quite straight- 
forward, so in the section devoted to the numerical 
implementation of the algorithm we focus on dis- 
cussing the ways to secure the algorithm robust- 
ness for the media with non-linear plasma density 
changing in space by presumably decades of or- 
ders of magnitude. This is shown to be the case 
for the heliosphere. We describe the methods we 
used to adaptively decrease and increase the inte- 
gration step and the techniques used to avoid the 
ray penetration into the regions with the density 
above the critical. These methods are intended 
to reverse the ray direction and use parabolic and 
linear approximations of the trajectory. 

In the testing and validation section we show 
the algorithm demonstrates excellent precision in 
the medium with one-dimensional linear plasma 
density gradient. In the example application sec- 
tion we verify the algorithm through solving a 
ray-tracing problem in the corona and chromo- 
sphere with realistic mod el density distributions , 
those bv lSaitol (|l97d) and lCillie fc Menzell (Il935l) . 
Our algorithm implementation shows good results, 
demonstrated in a few plots of the ray beams re- 
fracting near the sun. 



2. Fundamentals of ray refraction in isotropic 
plasmas 

We treat the electromagnetic ray trajectory 
as a curve with radius vector r in the three- 
dimensional space. For any starting point, rg, and 
any given vector of initial direction, vq, the whole 
shape of the ray trajectory may be unamnigously 
determined, as long as the distribution of the index 
of refraction, n{ r, lo), is known for a given angular 
frequency w. 

For the natural parameterization of the curve, 
r(s), where s is its arc length, the direction vec- 
tor v(s) is defined as the first derivative of po- 
sition vector, or v(s) = r'(s)- The components 
of v(s) are direction cosines, and |v(s)| = 1 at 
any point of the curve. Here we show that the 
plasma density distribution, p{ r), suffices to deter- 
mine the index of refraction n( r) at a given wave 
frequency. In the Gaussian units, n is the square 
root of plasma dielectric permittivity, n — ^/e. For 
isotropic plasmas 



n = 1 



(1) 



where the plasma frequency is = (47re^ne/me)"'^/^, 
e being the electron charge, TOe and the electron 
mass and number density in the Gaussian (CGS) 
units. Under the quasineutrality condition the hy- 
drogen plasma density is p = mpUe, where irip is 
the proton mass, so the dielectric permittivity is 
expressed as 

e = l-^ (2) 

Per 

Here pcr is the critical plasma density, at which 
the permittivity and, hence, the refraction index 
for the given frequency are zero: 



Pc 



47re2 



(3) 



Electromagnetic rays of the frequency u cannot 
travel in the regions with ujp > w, i.e. with p > 
Per where the index of refraction takes imaginary 
values. The boundary between the region where 
p < Per and the region with p > pcr is called the 
critical surface. The rays near the critical surfaces 
undergo strong refraction bending them away from 
the surface. In limit cases this effect is similar to 
the reflection. 
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The algorithm described further uses the rela- 
tive gradient of refractive index, Vn/n. Since the 
permittivity gradient can be expressed from Eq. 
(HI) as Ve = — Vp/pcr, the relative gradient of n is 



Vn 
n 



Ve 
"2e 



Vp 



2(pcr - P) 



(4) 



Thus, if we know the critical density, the relative 
gradient of refractive index is totally determined 
by the plasma density and the density gradient. 

3. Ray equation for isotropic media 

To derive the equation for the ray trajectory 
we use the Fermat's principle. It states that a ray 
joining two points Pi and P2 will choose the path 
over which the integral of the refraction index, n, 
reaches its minimum, so that the variation of the 
integral about this path is zero: 



n{ r)ds 



(5) 



Pi 



Notice that the integration in ([5]) occurs along a 
path, while in the Hamiltonian principle of least 
action the integration is with respect to the time. 
Transform the integrand of ([5]) into the form of 
the Lagrangian by changing s for a time-like inde- 
pendent variable r. If s = s(r), then ds^ = f^dr^, 
where = r • r. Substitution of As into (O, 



r-P2 

I n{ r)\/f2 dr = 0, 
Jpi 



(6) 



turns the integrand into the Lagrangian L{ r, r, r) = 
n(r)V7^. We apply the Euler method to ([6|) to 
find the vector differential equation of the curve 
satisfying the Fermat's principle. The Euler's 
equation is 



_d_ 

d7 



dL 
dr 



dL 
d r 



= 



(7) 



The Lagrangian derivatives with respect to r and 
r are 



dL , . r 
-— — n{r 



dr 
and 



/ X rdr , X dr 
^^(1-)^=^ =^(r)-T- (8) 



2 ^ ' ^/^2 ^ ' ds 



dL dn(r) v7 / \ 

—— ~ — vr^ = Vn(r)— — . 

d r d r dr 



(9) 



Therefore, the Euler's equation ([7]) for our prob- 
lem ([5]) takes the form 



d 

d7 



dr 
"dl 



Vn 



dr 



0, 



(10) 



which, after multiplying both parts by dr/ds, 
yields the ray equation: 



d. 



(11) 



With the use of the vector identities, Eq. pT|) can 
be represented in the format with explicit first and 
second derivatives of r as 

d^r '^'^ X X '^''^ (12) 

ds^ ds \ n ds J 

This is a system of three second-order equations. 
By introducing the direction vector v = dr/ds 
the system is transformed into the set of six 
first-order equations 




V X 



Vn 



X V 



(13) 



The set of equations ([T3| is numerically solved by 
the proposed ray-tracing algorithm. Note that the 
independent variable s is the ray arc length so the 
solution to (IT51) is a naturally parameterized curve 
r = r(s). This infers the "conservation law" for 
v(s) = r'(s), 

|v|^l, (14) 

at any point of the curve. The algorithm of ray 
tracing described in Section[S]is based on a scheme 
that conserves this property. 

4. Ray equation as particular case of 
Haselgrove equations 

It is intrersting to observe that equation ([TT|) 
is a particular case of a mor e gene ral set of ray 
equations derived bv lWalkerl (l2008l) : 



dr „ , , X 
— = Vfcw(r, k) 

G{uj, k) =0 



(15) 
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In Eqs. (fT5)) G (w, k) = is the dispersion relation 
of the medium, t is time, dr/df = r = Vg is 
the group velocity of the wave packet, uj is the 
angular frequency, k is the wave vector, and Vfc 
is the differentiation operator with respect to the 
k compone nts. Syste m (|15l) is a modification of 
the lHaselgrove I l|l955h equations. Its solutions are 
trajectories of rays of virtually any nature, not just 
electromagnetic, as long as the waves are plane 
and linear, varying in space and time as exp[i(k • 
r-ojt)]. 

Let us show that equation is equivalent 
to system (fT5|) for the waves traveling through a 
coUisionless isotropic plasma with the dispersion 
relation 

w2^c2/^2-Wp^ (16) 

which is easily found from equation dH) and the 
relation between the index of refraction and the 
magnitude k of the wave vector. 



n = ckjijj. 



(17) 



Hereafter c denotes the speed of light in the vac- 
uum. In equations ([T5|) we change variables, time 
t for arc length s. For a curve r = r(t), where t 
is a parameter, the element of arc length 



ds = \r\At = WcAAt. 



(18) 



Differentiation of ([T6| with respect to k gives the 
first equation of system ([TSj) : 



dr , 

— = Vfco; = — k, 
dr UJ 



(19) 



and therefore the group velocity and its absolute 
value are 

2 2 

Vg - Vfet^ = - k, Vg = -k. (20) 
Then, according to p^ . di is rendered as follows: 



(21) 



(22) 



The second equation of ()15p is obtained in a sim- 
ilar way by applying V to the dispersion relation 

(nni), 

dk _ uj„ 
"At 



At ~ -IT- ds. 

&k 

Substitution of At in equation ((T9)) yields 

dr k 

ds k 



-Vuj = — ^Vw, 
c 



(23) 



and using (|2T|) to change dt for ds: 

dk UJr, 



ds 



(24) 



We reduce equation ([22)) to a form similar to that 
of (|TT|) . Multiply both sides by ck/uj and use (|17l) : 



dr c , 
n— — = — k. 

ds UJ 

Differentiate both sides with respect to ds: 

d / dr\ c dk 

oj ds 



(25) 



ds V ds 



(26) 

Substitution of (p4)) in the right-hand-side yields 
d / dr\ „ ^27) 



ds V ds 



ujck 



The left-hand-side of (P7|) is the same as the first 
term in pT|) . while its right-hand-side is Vn. This 
can be verified by applying V operator to equation 
(U) and using the dispersion relation (IT51) . Thus 
we have shown that equation pip derived directly 
from the Fermat's principle in the previous Section 
is equivalent to the general ray equations (jl5|) con- 
strained by the conditions of plasma isotropy and 
absence of losses. 

5. Algorithm description 

The integration of system (IT51) is based on 
an original algorithm developed as a modifica- 
tion of the Boris' i mplem entation of the explicit 
Crank fc NicolsonI (jl947r ) scheme for the calcula- 
tion of a charged particle mo tion in magne ti c field 
(|Hocknev fc EastwoodI [l988h . The iBorisI (|l970l ) 
CYLRAD algorithm secures the energy conserva- 
tion, which is equivalent to conservation of the 
particle's squared velocity v^. The ray trajec- 
tory calculations require a different (though simi- 
lar) conservation law to be observed. 

Below is shown a brief derivation of the second- 
order difference scheme for the system p^ . We 
introduce a vector quantity ft as 

fl^ — xv (28) 
n 

Substitute p8|) into (|13p to render the ray equa- 
tions in the form 

d r 

"dJ ^ 
dv 

"d^ 



< 



(29) 



V X J7. 
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Formally, this set o f equations is the sam e as (4 - 
90) in the book bv iHocknev fc EastwoodI (|l988l ). 
However, in the ray trajectory equations (f29| the 
independent variable s, the ray arc length, replaces 
the time i, introducing the natural ray parameter- 
ization. The ray direction v(s) — r'(s) is used 
instead of the particle velocity v(t) = r'{t), and 
the vector n{s) is used instead of the cyclotron 
frequency fl — q'B/m . Analogously, the square 
of direction vector v^(s) = 1 must also be con- 
served in the course of numerical integration, be- 
cause v^(s) = = 1 for any regular point of 
a naturally parameterized curve. 

It is easy to show that is a vector whose 
length equals the ray curvature, defined as k = 
|r"(s)| for naturally parameterized curves. Ren- 
der the second equation in (f29|) as r" = r' x Jl. 
Multiply both sides of it by r'x on the left: 



r' X r" = r' X (r' x O) 



(30) 



Expanding the right-hand-side with the use of the 
vector identity and keeping in mind that r' • r' = 1 
and r' _L Jl, we get 



The second differential equation in (|29l) can be ap- 
proximated by the differences as 



vi - vq 

ds 



Vl/2 X ^1/2 



(33) 



Approxi mation of v^ /■?. as ^( vi 4- vq) yields the 
implicit ICrank fc NicolsonI ( 1947 ) scheme 



vi - vo = (vi + vo) X r2i/2 — (34) 

This scheme has the benefit of unconditional sta- 
bility. Multiply both sides of (p4l) scalarly by r2i/2 
to get the property: 



(35) 



Vi • fli/2 — Vq • r2i/2. 



which will be used later. To convert scheme (p4)) 
into the explicit one, express vi from ([34]) as 

ds 

Vl = Vq + (Vi + Vo) X r2i/2 — (36) 



and substitute it in the right-hand-side of ([34)) for 
itself. It yields 



n = -r' 



X r 



(31) 



But for the natural parameterization r' _L r", 
hence, I r2| — I r"|, and 



(32) 



The numerical solution is sought at the nodes of 
one-dimensional grid marked on the ray trajectory. 
Consider a grid whose nodes with respect to some 
i-th position are numbered as ... (i — 1), (i — 1/2), 
(i), -1-1/2), (i + l),... and so forth. We subscript 
a variable with if its value is taken at i-th point, 
with (1/2) for {i + l/2)-th point, and with 1 for 
{i + l)-th point. A leapfrog difference scheme is 
used, where the half-step points are used in cal- 
culations. Suppose a ray at a particular moment 
has its end position Tq and direction vq. The 
ray tracing scheme provides formulae to calculate 
its successive position ri and direction Vi after 
the step ds along the ray arc length. The vari- 
ables subscripted with 1/2 are intermediate and 
they are related to the middle of the step. Note 
that the density and its gradient are provided only 
once for each step at the step middle point so the 
(Vn/n) ratio is ascribed to the grid point {i + l/ 2). 



Vl = Vo + 2vo • Oi/2 — 



+ ^(vi + vo) X Oi/2 X rii/ 



2^ (37) 

Using the identity Vo = vo(l -I- — 
{ill/2 -x)^)j vector identities, property ([55]) . and 
rearranging the terms, we arrive at the explicit 
difference scheme know n as the Boris' CYLRAD 
algorithm (Boris 1970[) : 



Vl 



Vo 



l + (ni/2f)' 



_ ds\ _ ds , , 
Vo -I- Vo X ni/2 — X ni/2 — (38) 



The vector r2i/2 requires extrapolation: as seen 
from (f28| . the first multiplier, Vn/n, is obtained 
at the intermediate point and therefore has the 
(1/2) subscript, while the second, v, is only avail- 
able for the starting point subscripted with 0. We 
denote this cross-product as fto, 



rio = I — I X Vo, 



(39) 
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and use it as an initial value for $7 extrapolation 
up to r2i/2- The half-step difference scheme for 
the second equation in ([29|) is 



Vi/2 = Vo + Vo X 



ds 



(40) 



Substitution of 'Vi/2 into (|28l) yields the extrapo- 
lated r2: 



^1/2 — ( 



Vo X ^0 ^ ) (41) 



In order to make the algorithm faster we reduce 
the number of multiplications. Note that all the 
occurrences of ft in ((38|) and (|4T|) have the mul- 
tiplier ds/2, so it is expedient to include it in the 
expressions for in (|39p and (|¥T]) . Finally, we put 
together the formulae (jSH) , (HTI) , and (|55)) and add 
the position vector r adjustments at the start and 
the end to form the ray tracing algorithm (|42|) . 

The algorithm presented below converts the ini- 
tial ray position and direction (ro, vq) into the 
new ones ( ri, vi) along the ray path in five steps: 



ds 

1. ri/2 = ro+ v„ — 

„ ^ /Vn\ ds 

2. r2o = X Vo — 

\ n I 1 2 



3. rji 



Vn 



/2 



4. Vi = Vq 



X (vq + Vo X f2o) „ 

nil 2 



(42a) 
(42b) 
ds 



(42c) 



'1/2 — 

(vo vo X ^1/2) X f2i/2 (42d) 
ds 

5. ri = ri/2 + vi — (42e) 



The first step (|42ap provides the r half-step- 
ahead approximation Next two steps (|42bp 
and (|42cp calculate the half-step-ahead approxi- 
mation for SI. The fourth step (j42dp estimates the 
new ray direction, vi, based on the previous di- 
rection, Vo, and ^1/2- Step five, (|42el) . produces 
the new ray position, ri, through the half-step 
from the intermediate point Vi/2 along the new 
direction vi. 



6. Numerical implementation 

The optical properties of the media where trac- 
ing the electromagnetic rays is of interest are often 
extremely nonuniform. For example, the electron 
number density Nf, within the chromosphere alone 
falls off with the height by more than 10 orders of 
magnitude. Inside the solar corona decays with 
the solar distance also by decades of magnitude or- 
ders. The index of refraction, as the square root 
of proportional to N,, dielectric permittivity ([2]), 
can change by many orders of magnitude within 
the region of study. The non-linear dependence of 
the refractive index on the wavelength poses other 
issues. For different wavelengths the ray geome- 
try can be very different because the regions with 
the density greater than critical ([3]) form differ- 
ent configurations in space. The numerical imple- 
mentation of algorithm therefore must include the 
mechanisms for adaptive step size ds adjustment 
not only to ensure the specified precision but also 
to prevent the ray from penetration into the re- 
gion with p > Per- In some cases ro appears so 
close to the region with the critical density that 
the step ds must be shorter than the absolute min- 
imum step specified. For this condition a special 
method of replacing the ray trajectory by a seg- 
ment of parabola or even straight line is used. All 
these techniques are described in this Section. 

In our implementation of the ray tracing al- 
gorithm, aimed at obtaining multi-pixel images, 
many rays are processed simultaneously. Over a 
single integration step each of the rays is tested 
for precision and correctness, and only those that 
passed the tests are submitted to the Boris's algo- 
rithm to calculate new positions and directions. 
The rays for which the absolute minimum step 
leads to penetration to the region with p > pcr 
are switched to the opposite parabola branches or 
just reflected from the critical surface. The rest 
are the rays whose step sizes are too large to se- 
cure the specified precision. The step size for them 
is shortened by applying formula (|46p . The rays 
with reduced step sizes are processed at the next 
integration step. If a ray is leaving the dense re- 
gion, its step size ds is increased over the series of 
subsequent integration steps. 

The plasma density, p, and its gradient, Vp, 
are calculated through calls to the specialized sub- 
routine once per each algorithm step. It is as- 
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sumed that these plasma parameters are obtained 
through an MHD simulation and they are avail- 
able at the nodes of a three-dimensional grid with 
variable grid step. The same subroutine provides 
the estimate for the maximum ray step size con- 
sistent with the MHD grid size for each ray. 

6.1. Precision control through adaptive 
step size 

The desired precision of the ray trajectory is 
specified with the variable Tol . It has the meaning 
of inverse number of the steps Us per a radian of 
the ray arc length. For example, Tol < 0.1 means 
that at least 10 steps per one radian of ray arc 
length is required. At each algorithmic step the 
real precision is calculated and compared with Tol, 
and if it exceeds Tol, the step ds is decreased. 

The curvature k is defined as the inverse of in- 
stant curvature radius as k = 1/Rc- The actual 
ray curvature is available at each algorithm step 
as I n| due to ([5^ . Instead of the real precision 
it is convenient to use the quantity r^, which is 
actually a halved real precision: 



to (1431) it is corrected as: 



n 



ds 



ds 



(43) 



The algorithm step (I42dp requires calculation of 
anyway, so using it saves time. It is that is 
compared with ToP to make decision on the ds re- 
duction. Comparing the halved real precision with 
the desired precision Tol makes the step ajustment 
less frequent as the specified precision is secured 
over most of the ray trajectory portions. To re- 
duce the number of floating-point calculations, the 
approximation uses the condition = 1 in the 
identity ( ux v)^ = u^v^ — (u-v)^= u^ — (u-v)^. 
The index of refraction n in ([^5]) is replaced by the 
dielectric permittivity as defined in (|4]), so 



ds /Vei/2 



X vo 



2 V 2ei/2 

which with the use of the identity yields 



(44) 



ds 
T 



/2 



Vo 



2ei/2 / \ 2ei/2 

(45) 

The step must be ds — Rc/rig ~ ToI/k to have 
exactly I /Tol points per radian of the ray, so due 



ds' 



Tol ds 



2 ' 



(46) 



The subroutine that provides the plasma den- 
sity and density gradient at each call also provides 
the desired new step size ds„eto for each ray. Af- 
ter the positions and directions for the rays have 
been advanced with the current step sizes, the al- 
gorithm attempts to increase the steps towards 
dsnew However, ds„eu; can be greater than the 
current ds by orders of magnitude, and near the 
critical surface a simple assignment ds — dsnew 
would almost inevitably lead to the ray penetra- 
tion past the critical surface. Therefore the algo- 
rithm uses a special method for gradual increase 
of the step size. It is based on running a non-linear 
difference scheme which outputs the increased ds' 
from the current ds and the desired dsmw 



Solutions to y^^^={2-yjx^)y. 




Fig. 1. — Solutions to the non-linear difference equa- 
tion H47p for four different initial conditions yo. The 
input variable is constant: Xi = 1. The solutions, yi, 
are increasing functions, always converging to Xi = 1. 
This scheme is used to smoothly increase the step ds 
to the set value of ds,ie™. 

Consider a difference equation associating two 
variables, the input x and the output y, specified 
at the points ... (i-1), (i), (i+1), (i+2) ... and so 
on: 

Vt+i = (2 - —)yi 



(47) 



For constant input Xi = X the solution to Eq. (|T7)) 
is stable for any initial condition on the output 
variable t/o between and X. The solution to 
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Eq. (|47)) always converges to the value X. The 
useful property of Eq. (|47|) is that for a small 
(compared to X) initial condition yo the solution 
behaves similar to growing exponential function 
e*, while for large values of yo, i.e. those closer 
to the X, it grows slower and slower similar to the 
decaying exponent Xe~*, eventually converging to 
X. This behavior is illustrated in Fig. [T] Notice 
that reasonably large initial values for yo ensure 
fast convergence within 1% of the range [0..1]: 10 
steps for yo = O.OIX, and even 3 steps for yo = 
0.5X. 

In our algorithm implementation the step ds is 
incremented using the Eq. (|47p as the algorithm 



ds' = 2 - 



ds 



(48) 



When ds is small compared to dsnew, each 
next value of ds' is almost 2ds, so ds grows 
in a geometrical progression. However, as ds 
approaches dsnetu, its growth rate becomes slower. 
These features ensure fast and yet smooth adjust- 
ment of ds in several steps up to the specified 
ds„ei« value. 

6.2. Correctness control 

Correctness of the ray path computation is vi- 
olated if at some step the ray penetrates the criti- 
cal surface. This condition can arise for very steep 
rays, almost normally incident on the critical sur- 
face. To prevent a ray from trespassing the critical 
surface the algorithm makes a linear prediction on 
the next value of dielectric permittivity e. If the 
predicted e is negative, the ray is approximated by 
a parabola and it is then switched to the symmet- 
rical point on the opposite branch of the parabola. 
This is similar to the reflection at the critical sur- 
face, after which the ray starts traveling away from 
it. However, the linear e prediction may fail due to 
highly non-linear plasma density distribution and 
then the ray trespasses the critical surface. This 
unphysical condition is corrected by returning the 
ray its previous position and direction, numerical 
calculation of the distance to critical surface along 
this direction, and linear reflection of the ray at 
the critical surface. 

When the current ray state is (tq, Vq) and the 
algorithm is calculating its next position and di- 
rection ( ri, Vi), the mid-step values of the plasma 



density pi/2 and its gradient Vpi/2 are obtained 
from the subroutine, so the dielectric permittiv- 
ity £1/2 > and its gradient Vei/2 are available 
from ^ and respectively. If the ray at next 
step occasionally trespasses the critical surface, 
the next value of permittivity, £3/2, will be neg- 
ative. Therefore, predicting the sign of £3/2 helps 
the algorithm to make decision whether to use the 
Boris' scheme (j42|) or take measures to avoid the 
penetration past the critical surface. The predic- 
tion is based on testing the condition 



£3/2/£l/2 > 



(49) 



Replace both numerator and denominator in 
49} by the linear portions of their Taylor series: 



£3/2 _ £0 + Vo • V£i/23 



ds 



£1/2 £0 + Vo • V£i/2 ^ 
Vo • V£i/2 



= 1 



£0 + Vo • V£i/2 ^ 



ds > 0, (50) 



where £0 is the permittivity value at the point 
ro . After rearranging the terms and dividing both 
parts by £0 we arrive at the correctness criterion: 



ds V£i/2 \ ^ 1 

" £0 



(51) 



Since the plasma density subroutine provides only 
the half-step value, pi/2 and, hence, only the half- 
step permittivity value, £1/2, can be directly cal- 
culated, £0 is obtained as its back approximation 
using the linear terms of its Taylor series: 



£0 



£1/2 - Vo • V£i/2 



ds 



(52) 



If inequality (|5T|) is true, the next step is supposed 
to be safe. If not, the ray is so close to the critical 
surface and so steep that it will penetrate the crit- 
ical surface at the current step ds. In this case the 
next position and direction are calculated through 
parabolic approximation. 

6.2.1. Parabola switching 

This method is used when the correctness crite- 
rion (|5T|) fails and the ray is about to trespass the 
critical surface. If the step ds is small enough, the 
dielectric permittivity gradient V£ can be consid- 
ered constant over the ray path. In the medium 
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with Ve — const the electromagnetic rays have 
parabohc trajectories, which will be shown fur- 
ther. Knowing the parameters of the parabola 
(see Fig. [5]), the algorithm calculates the posi- 
tion and direction vectors (ri, Vi) of the ray at 
the symmetric point of the opposite branch of the 
parabola. The ray is switched to the opposite 
parabolic trajectory branch so that instead of ap- 
proaching the critical surface it starts departing 
it. In order to find the new position and direction 
after the parabola switching we solve the ray Eq. 
([TT|) in two dimensions in the XF-plane contain- 
ing the vectors — Veo and vo, the X axis collinear 
with — Veo- 




L X 

Fig. 2. — In a medium with linearly changing dielec- 
tric permittivity e electromagnetic ray trajectories are 
second-order parabolae. The e linearly changes from 
1 at a: = to at a: = L. The critical surface is the 
locus where £ = 0. 

Consider a simple 2-dimensional case of a ray 
traveling within the XK-plane, where the dielec- 
tric permittivity changes linearly from e = 1 at 
the origin of coordinates down to e = at some 
distance L along the X axis, so that 

£ = l"f (53) 

The case schematic is shown in Fig. [2] We shall 
derive the equation for the ray trajectory passing 
through an arbitrary point (a;o,?/o)- The equa- 
tion set for the ray trajectory ([TT|) in the two- 
dimensional geometry consists of two equations for 
X and y. The dielectric permittivity e changes in 
the X direction only, so does the index of refrac- 
tion n, therefore its gradient component along the 



Y axis, dn/dy, is zero, and the equation for y takes 
the form 

which immediately implies 

n — — — const (55) 
as 

Refer to the schematic in Fig. [S] Denote the values 
of n and e at the point (xq, yo) as no — where 
£o is calculated as in ([52l) . The direction vector 
component dy/ ds — sin a, where a is the angle 
between the ray direction vq and — Veoj which is 
the X axis here. Therefore, the constant in (|55|) 
can be presumed in the form uq sin a, so 



Using the relationship dx^ / ds"^ + dy^/ds^ — 1 
and excluding ds, we arrive at the differential 
equation 




Fig. 3. — Parameters of a ray trajectory approxima- 
tion by the quadratic parabola in the XY plane. The 
X axis is aligned with the negated dielectric permit- 
tivity gradient, — Ve. The direction cosines' vector vo 
at the point (a;o,j/o) is decomposed into two compo- 
nents, V|| along the X and vx along the Y, with the 
lengths equal to cos a and sin a, respectively. When 
the ray leaves the parabola at the point {xi,yi), its 
X coordinate remains the same, but y changes by Ay, 
equal to the doubled distance from yo to the parabola 
vertex y* along the Y axis. 
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Since e is a linear function ([53]) . e{x) = Eq + 
{x - xo)£o, where e'q = de/ dx\x=xo, so Eq. ([57]) 
can be integrated as 

dx dy 



y/ EqX + EqXo - eo cos^ a Jy^ V^osma 

(58) 

to yield the solution in the form x = x{y): 



^0 2 
X = Xq y cos Q!+ 



1 / ^ e; 

— VEocosa 



" (y-2/o)V (59) 



2-v/eo sin a 



This second order parabola equation describes 
the electromagnetic ray trajectory for linearly dis- 
tributed e. Its vertex point {x*,y*) coordinates 
are 

y* = yo — —r cos a sin a (60a) 

^0 

X* = xq — cos^ a (60b) 





Vi = 


Vo - 2v|| 






■-yi 



Fig. 4. — The ray direction vi at the opposite branch 
of the parabola after the parabola switching or Snell's 
reflection. 

We intend to avoid tracing the points lying on 
the parabola. Assume the ray has already traveled 
from the point (xo,?/o) Bill the way to the axially 
symmetric point (a;i,yi) at the opposite branch. 
In three-dimensional case its position and direc- 
tion have changed to (ri, Vi). We need to find 
the increment A r such that 



ri ro + Ar 



(61) 



Note that |Ar| = Ay = 2{y* - yo) and = 
— jVeol because e' < in the considered vicinity. 
Using (|60al) yields the length of the increment: 



lArl 



4£o 



■ cos a sm a 



(62) 



In the chosen XY coordinates A r points in the Y 
direction. Decompose vq into its X and Y com- 
ponents V|| and v±, 



Vo 



^11 



(63) 



Here | Vo| = 1, | V|| | = cosa, and | vj_| = sina. We 
can state that the vector increment A r should be 
made in the Y, i.e. v_l, direction, so 



Ar = lArl 



Since a is the angle between — Veo and vq, 

Vo • (-Veo) 



nil 



iVeol 



(64) 



(65) 



whence the vector increment from ro to ri at the 
opposite branch of the parabola is 



A A I ^0- Veo ' 
A r = 4e I — — — — \ vj_ 



(66) 



From (|53p the (hypothetical) distance to the 
critical surface is L = l/|Veo|. With cosa from 
()65p . we introduce a useful expression 



L cos a 



Vo • Veo 



(67) 



iVsoP ' 

using which the increment (j66|) can be written as 

A r = 4eL cos a vj^ (68) 

The new ray direction Vi is calculated from the 
geometrical considerations shown in Fig. 2] as 



Vi = vq — 2 vi 



(69) 



To obtain the vector V||, we multiply its abso- 
lute value I V|| I from (|65p by the unity vector along 
the X axis, which is — V£o/|Veo|, to get the result: 



VfW^Veo) 
2 — ^^0, 



|Veo 

or, using the (p7)) notation, 

V|| = — LcosaVeo 



(70) 



(71) 



Thus, after the parabola switching, the new ray 
position and direction are 

ri = ro + ieL cos Q!( vo + L cos aVeo) (72a) 
vi = Vq + 2L cos aVeo ('''2b) 
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{xi,yi) 








Ax 


(xo, yo) 



Fig. 5. — The length of parabohc segment ABC is 
approximated by the added lengths of two parabolic 
chords AB and BC. This length equals that of the 
line segment AD. 

For the purposes of path integration in line with 
the trajectory calculation, the integrand is multi- 
plied by the length of every step ds. In case of 
the parabola switching the integrand is multiplied 
by the length of the parabola segment, dsp, from 
ro to ri. Assuming this length is small enough 
it can be replaced by the length of two straight 
line segments connecting the points (in XY coor- 
dinates) (a;o,2/o): the parabola vertex (a;*,y*), and 
{x\,y\), as shown in Fig. [5l Denote as Aa; the dis- 
tance from (zq, yo) to (a;*, y*) along the X axis or, 
which is the same, along the vector — Veq (see also 
Fig. [3]). Then the parabola length approximation 

is 

dsp = V(2Aa;)2 + Ay^ (73) 

We already know that Ay — |Ar|. From (j60bp 
and using (|67| we can find: 

Ax = e(Lcosa)^|V£o|- (74) 

6.2.2. Linear reflection 

The criterion (I5l1) is based on the linear for- 
ward approximation of e, and sometimes it works 
incorrectly if the plasma density distribution near 
the critical surface is highly nonlinear (for exam- 
ple, exponential, as in the chromosphere). Under 
these conditions the calculated new point of the 
ray, ri, appears in the region with p > pcr, which 
is impossible in the physical reality. To fix his con- 
dition the algorithm stores in memory the previous 
ray position and direction (r_i, v_i). 

When the penetration is detected, the algo- 
rithm restores the ray state to one step back: 

(ro, vo) 4- (r_i, v_i). (75) 



The parabola switching mechanism in this sit- 
uation can not be relied upon because the lin- 
ear Taylor series approximations similar to ()52|) 
have already failed and the parabola parameters 
would be calculated incorrectly. The algorithm 
uses the Newton's method to find the length of 
step ds along the direction Vq such that the next 
ray point, ri, lies on the critical surface with 
some small clearance Tolg. Then the algorithm 
performs linear reflection according to the Snell 
law. Since this conditions only occur at high fre- 
quencies when the rays are close to straight lines, 
this procedure does not deteriorate the algorithm's 
precision. 

The Newton's algorithm is based on stepwise 
approximation: 

dsi+i = dsi - ■=r^: (76) 

where Vue^ = Vq ■ Vci is the directional derivative 
of Ei along the vector Vo, so the calculations are 
performed according to the scheme 

dsi+i = dsi (77) 

Vo • Vci 

The process stops when |ei| ^ Tol^. The ray di- 
rection is then corrected according to the Snell's 
law. Its geometry is the same as for the parabola 
switching shown in Fig. |4l so 

Vi=vo-2v|| (78) 

where V|| is calculated as in ([7T|) . and the next al- 
gorithmic step is performed with the last ds value 
immediately from the critical surface. 

7. Testing and validation 

The software implementation of the ray trac- 
ing algorithm comprises computer programs in the 
Fortran 90 and C languages. As an example of 
the algorithm validation, the ray trajectories in 
a hypothetical medium with the linear dielectric 
permittivity e distribution have been obtained nu- 
merically using the algorithm and compared with 
the exact solutions. The results are shown in Fig. 

n 

The permittivity e linearly falls off along the 
X axis starting from 1 at x = down to at 
some distance i, as given by Eq. (j53p . The ray 
trajectories in this medium have been shown to 
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Fig. 6. — Comparison of the numerically simu- 
lated ray trajectories for several angles of incidence, 
l^jlO^jSO", and 60'', with the exact ray traces. In 
panel a the simulated ray curves, shown in thick 
dashed lines, overlap the exact parabolic solutions, 
plotted in solid thin lines. Panel b shows a closer 
view of short segments of the 30° traces near the re- 
flection point at a; = 75 calculated by the ray tracing 
algorithm (segmented line) and the exact solution H79[) 
(smooth, solid line). Notice the scales of the axes. 



be parabolas described by Eq. ([5^ . We chose the 
distance L — 100 of conditional units (they might 
be solar radii, for example), and set the starting 
point To at the coordinate origin (0,0,0). Here 
the ray equation (|59p takes the form 



X — L cos a 



1 



cos a 



2L sin a 



(79) 



where a is the angle of incidence on the critical 
surface at x = L. The results of comparison are 
presented in panel a of Fig. [51 Four angles a were 
selected: 1°, 10°, 30°, and 60°. The exact solutions 
calculated by ([7^ are plotted in thin solid lines; 
the ray traces created by the algorithm are plot- 
ted over them as thicker dashed lines. The dashed 
curves cover the solid ones, so one can conclude 
the numerical solutions are essentially precise. In 
panel b of Fig. IH] the differences between the nu- 
merical and exact solutions are more discernible 
because only a short piece of the trajectory in the 
vicinity of the parabola vertex is shown. 

8. Example applications 

The algorithm has been applied for solar stud- 
ies and shown good results. A series of numerical 
experiments was conducted to simulate the bright- 



ness temperature distribution, Tg, across the solar 
disk as seen from the earth, in the wide frequency 
range from 10 MHz to 3 GHz. The calculation 
of Tb requires the ray tracing. We used different 
electron density distributions, Ne (cm~^), for the 
corona and the chrom osphe r e. Th e coronal Ng was 
approximated by the Saitol ( 1970f l model: 



Ne{r,9) =3.09x 10V"i°(l - 0.5 cos0)+ 
1.56 X 10V"°(1 - 0.95 cos 6*)-^ 



0.0251 X 10V"2-^(1 - Vcos0), (80) 

where r is distance in solar radii, Rq, from the 
sun center and 6 is the heliographic colatitude. In 
the chromosphere, up to 9,000 km from the photo- 
sphere, we used the lCillie fc Menzell (|l935l) model: 



iVe(r) = 5.7 X 10 



llg-7.7xl0" 



-500) 



(81) 



and in the layer between the chromosphere and 
corona, for the altitudes from 9,000 km to 11,000 
km, a smooth polynomial patch function was used. 

Electromagnetic Rays Refracting Near Sun at 80.0 MHz 




Fig. 7. — The result of numerical simulation of a beam 
of electromagnetic ray traces at 80 MHz in the solar 
corona. The square cross-section beam consists of 7 x 7 
evenly spaced rays direct ed along the X axis i n the 
HEE coordinate system (jPranz fc Harpeij|2002l ). i.e. 
from the earth. The sun is shown as a yellow sphere. 

In Fig. [7] is shown a bunch of the computed ray 
trajectories at 80 MHz, refracting in the coronal 
and chromospheric plasmas. The 7 x 7 = 49 evenly 
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spaced rays comprise a beam converging to a point 
at the earth distance (here 215 Rq), where the 
rays can be displayed as a 7 x 7-pixel raster image. 
Note that in our experiments with the ray-tracing 
algorithm the typical raster sizes were from 100 x 
100 to 1000 X 1000 and more; the beam in Fig. [7] 
is shown for illustrative purposes only. As seen, 
the rays arriving at the earth to form an image at 
80 MHz can carry information from very different 
heliospheric regions due to differences in refraction 
across the image plane. 

Fig. [8] shows the results of ray tracing at 
200 MHz in the XZ plane of the the HEE 
( Franz fc Harper 2002[ ) coordinate system (the X 
points from the sun center at the earth). The rays 
make up an envelope outlining the critical sur- 
face, which is at approximately 1.2_Rq from the 
sun center. 



Electromagnetic Rays Refracting Near Sun at 200.0 MHz 
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Fig. 8. — Refraction in the corona at 200 MHz 
in the XZ plane of the HEE coordinate system 
(|Franz fc Harped 1200 j ). The X points at the earth, 
the Z is normal to the ecliptic plane. The set of prop- 
agating rays is shown in thin solid lines. The zone of 
avoidance close to the solar surface, outlined by the 
envelope of rays, is clearly visible at 1.2J?0. 

The higher the frequency, the deeper the elec- 
tromagnetic rays penetrate into the solar vicinity. 
Fig. [9]helps estimate the solar distances where the 
rays at different frequencies start to significantly 
change the direction. The two thin beams of rays 
are directed parallel to the HEE X axis at the 
Z distances of O.18i?0 and O.67i?0, respectively. 
Each beam includes rays at the frequencies 10, 18, 
40, 80, 200 MHz, and 3 GHz. The Fig. Hcan be 
also viewed as the "spectral decomposition" . One 



can notice that at 3 GHz the rays are virtually 
straight lines. 



Refraction Near Sun for Several Frequencies 
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Fig. 9. — Refraction near the sun for several fre- 
quencies ("spectral decomposition"). Two thin beams, 
each composed of the frequencies 10, 18, 40, 80, 200 
MHz, and 3 GHz, are directed along the X axis in the 
XZ plane in the HEE coordinates l|Franz fc Harpej 
I2OO2I ). The X points at the earth, the Z is normal to 
the ecliptic plane. The beams' axes are at 2 = 0.18-Rq 
and z = O.67-R0. One can see the observations at dif- 
ferent frequencies are dominated by different layers of 
the corona. 



9. Conclusions 

The proposed algorithm is intended for efficient 
computation of electromagnetic ray traces in the 
heliospheric plasmas in parallel with integration 
over the ray paths. To be fast, it is based on a 
difference scheme of the second order of precision, 
but it has a property of conserving the direction 
cosine vector length of the computed trajectory. 
With this benefit, the algorithm maintains suf- 
ficient precision. The difference scheme is only 
a part of the algorithm: besides, it incorporates 
several mechanisms for adaptive estimation of the 
optimal step and prevention of occasional ray pen- 
etrations into the regions where the plasma fre- 
quencies are higher than the own ray frequency. 

The algorithm has been tested on problems of 
the solar disk image simulation, and it demon- 
strated good speed, precision, and robustness. We 
presented here some results pertaining to the ray 
traces calculation only; the main results will be 
published in a separate paper. However, we can 
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conclude that, as seen from Fig. |9l observations 
at different frequencies are dominated by different 
layers of the corona and chromosphere The broad 
band observations could help reveal the large scale 
structures in these domains, which can be likened 
to peeling the layers of an onion bulb. The new 
high-resolution radiotelescope MWA with its fre- 
quency range from 80 to 300 MHz therefore will 
be an effective instrument in the solar studies. 
The algorithm and software described in this pa- 
per were intended for use as the analysis tool for 
the MWA solar images. 

However we should specially point out that the 
algorithm applications are in no way confined to 
the MWA project only. The algorithm can be used 
for virtually any radioastronomical study because 
of its very wide frequency range. 

Among possible future improvements to the al- 
gorithm we could mention the plans on employ- 
ment of the multi-core, hyperthreading, and multi- 
channel-memory powers of the modern comput- 
ers. The problem of simultaneous computation 
of many ray trajectories can be effectively decom- 
posed into several threads, each running on a sep- 
arate CPU, because the rays in a beam are calcu- 
lated independently. 
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